################################################################################
###
### Filename: 	Maps-Slavery-African-Americans.R
### Author: 	  David A. Bateman
### Function: 	Reproduces maps of the distribution of enslaved African 
###             Americans and later of all African Americans shown in 
###             online appendix of Deeper Roots.
###             Parts of this code are based on slavery-maps-replication.R, 
###             submitted by Acharya, Blackwell, and Sen to the Journal of 
###             Politics dataverse to reproduce figures and tables in Acharya, 
###             Blackwell, and Sen (2016). 
###	Last update: January 1, 2023
###
################################################################################
rm(list = ls())

library(readxl)
library(dismo)
library(ggmap) 
library(raster)
library(rgdal)
library(rgeos)
library(elevatr)
library(plyr)
library(reshape)
library(ineq)
library(maps)
library(maptools)
library(foreign)
library(RColorBrewer)
library(classInt)
library(shapefiles)
library(ggplot2)
data(state.fips)
state.fips <- unique(state.fips[,c("fips","abb")])
state.fips$abb <- as.character(state.fips$abb)
state.fips <- rbind(state.fips, c(2, "AK"))
state.fips <- rbind(state.fips, c(15, "HI"))
rownames(state.fips) <- state.fips$abb
fips.state <- state.fips
rownames(fips.state) <- fips.state$fips
data(county.fips)

dir <-paste("../Deeper Roots/", sep="")
setwd(dir)

################################################################################
### Historical county data available from http://publications.newberry.org/ahcbp/
### Load county and state shape files. Identify the South using those states
### included in Acharya, Blackwell, Sen (2016).
histcounties <- readOGR("Dataverse/Shapefiles/US_AtlasHCB_Counties_Gen001/US_HistCounties_Gen001_Shapefile", layer = "US_HistCounties_Gen001")
histcounties$fips <- as.numeric(as.character(histcounties$FIPS))
histstates <- readOGR("Dataverse/Shapefiles/US_AtlasHCB_StateTerr_Gen0001/US_HistStateTerr_Gen0001_Shapefile", layer = "US_HistStateTerr_Gen0001")
ss2000 <- which(histstates$START_N <= 20001231 & histstates$END_N >= 20001231)
states2000 <- histstates[ss2000,]
states2000 <- states2000[!(states2000$ABBR_NAME %in% c("AK", "HI")),]

southern.states <- c("Alabama", "Arkansas", "Florida", "Georgia", "Kentucky", "Louisiana",
                     "Mississippi", "Missouri", "North Carolina", "South Carolina", "Tennessee", "Texas",
                     "Virginia", "West Virginia", "Maryland", "Delaware")


south <- states2000$ABBR_NAME %in% c("AL", "AR", "MS", "LA","TX","FL","MO","KY","WV","VA","MD","DE","NC","SC","GA","TN","DC")
states2000$south <- south
the.south <- unionSpatialPolygons(states2000, south)

## These data sources are originally from Jeremy Atack. The links and code are 
## from interpolation.R (Acharya, Blackwell, and Sen (2016)).
## https://my.vanderbilt.edu/jeremyatack/data-downloads/
rail <- readOGR("Dataverse/Shapefiles/RR1826-1911Modified0509161", layer ="RR1826-1911Modified050916")
rivers <- readOGR("Dataverse/Shapefiles/SteamboatNavigatedRiverContigUSAprojection", layer = "SteamboatNavigatedRiverContigUSAprojection")
canals <- readOGR("Dataverse/Shapefiles/19thC_Canals_March2017", layer = "19thC_Canals_March2017")

### Create heat map function.
plot.heat1 <- function(counties.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map, border=NA, lwd=1,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat), add=TRUE)
  
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}
color.palette1 =  colorRampPalette(c("white","grey","black"),space="rgb", bias=1)

### Make sure they're all on same projection.
counties.map <- spTransform(histcounties, CRS(proj4string(rivers)))
canals <- spTransform(canals, CRS(proj4string(rivers)))
the.south <- spTransform(the.south, CRS(proj4string(rivers)))

################################################################################
### Create maps for 1830, 1850, 1860, and (using African American population)
### for 1900. Instead of a loop, we just set the year by assigning it to 
### year1. 
################################################################################
year1 <- 1900
year2 <- year1+10 
year3 <- ((year1-1)*10000)+101

### The data need to be set by decade. The rail data doesn't have entries
### for 1830 so don't project if that's the case.
rail.x   <- rail[rail$InOpBy < year1,]
rivers.x <- rivers[rivers$START < year1 & rivers$END > year2,]
canals.x <- canals[which(canals$OPENED < year1 & canals$CLOSED > year1),]

if(year1==1830){
  fig<-12
}
if(year1==1850){
  fig<-13
  rail.x <- spTransform(rail.x, CRS(proj4string(the.south)))
}
if(year1==1860){
  fig<-14
  rail.x <- spTransform(rail.x, CRS(proj4string(the.south)))
}
if(year1==1900){
  fig<-15
  rail.x <- spTransform(rail.x, CRS(proj4string(the.south)))
}

## Data comes from ICPSR 35206
## http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies/35206
file.x <- paste("Dataverse/Census/census",year1,".dta",  sep = "")
census.x <- read.dta(file.x)
vartable <- cbind(names(census.x), attr(census.x, "var.labels"))
if(year1<1862) {
  cvars <- c("totpop", "stot", "fips")
  census.x <- census.x[census.x$level == 1, ]
  census.x <- census.x[, cvars]
  census.x <- census.x[!is.na(census.x$fips),]
  census.x$fips <- floor(census.x$fips)
  census.x <- cast(melt(census.x, id = "fips"), fips ~ variable, sum)
  census.x$xprop <- 100*(census.x$stot/census.x$totpop)
}
if(year1>=1865) {
  cvars <- c("totpop", "colmtot", "colftot", "fips")
  census.x <- census.x[census.x$level == 1, ]
  census.x <- census.x[, cvars]
  census.x <- census.x[!is.na(census.x$fips),]
  census.x$fips <- floor(census.x$fips)
  census.x <- cast(melt(census.x, id = "fips"), fips ~ variable, sum)
  census.x$xprop <- 100*((census.x$colmtot + census.x$colftot)/census.x$totpop)
}

#counties.map <- spTransform(histcounties, CRS(proj4string(rivers.x)))
#canals.x <- spTransform(canals.x, CRS(proj4string(rivers.x)))
#rail.x <- spTransform(rail.x, CRS(proj4string(rivers.x)))
#the.south <- spTransform(the.south, CRS(proj4string(rivers.x)))

### Some of the counties need correcting (especially SC, which used districts)
### Most of this code is from Acharya, Blackwell, Sen (2016)'s interpolation.R.
### New Mexico also poses problems in 1850 - drop
counties.map$fips[which(counties.map$NAME == "CARROLL (ext)")] <- 22035
counties.map<-counties.map[counties.map$STATE_TERR != "Mexican Cession",]
counties.map<-counties.map[counties.map$STATE_TERR != "New Mexico Territory",]
cc.x <- which(counties.map$START_N <= year3 & counties.map$END_N >= year3)
counties.x <- counties.map[cc.x,]
kk <- as.character(counties.x$NAME[counties.x$STATE_TERR == "South Carolina"])
kk <- sub(" Dist.", "", kk)
kk <- paste("south carolina,",tolower(kk), sep = "")
kk.fips <- NA
kk.fips[kk %in% county.fips$polyname] <- unlist(lapply(kk, function(x) county.fips$fips[county.fips$polyname == x]))
counties.x$fips[counties.x$STATE_TERR == "South Carolina"] <- kk.fips
counties.x <- counties.x[!is.na(counties.x$fips),]
newpolys <- unionSpatialPolygons(counties.x, counties.x$fips)
kk <- unique(as.data.frame(counties.x)[,c("fips", "STATE_TERR")])
rownames(kk) <- as.character(kk$fips)
counties.x <- SpatialPolygonsDataFrame(newpolys, data = kk)

counties.x@data$order <- 1:nrow(counties.x@data)
counties.x@data <- merge(counties.x@data,census.x, by=c("fips"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.x@data<-counties.x@data[order(counties.x@data$order),]


################################################################################
### Figure A12-A15: Distribution of enslavement/African Americans            ###
###             in the South, across different years                         ###
################################################################################
label<-paste("Dataverse/Figures/Figure-A", fig, ".pdf", sep="")
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(1,0,1,0), mar=c(0,0,2,0))

plot(the.south[2,], density=c(10))
plot.heat1(counties.x[counties.x$STATE_TERR %in% southern.states,],z="xprop",breaks=c(0:100), plot.legend=FALSE)

plot(the.south[2,], add=TRUE)
plot(rail.x, col = "indianred", add = TRUE)
plot(rivers.x[rivers.x$END > year1,], col = "dodgerblue", lwd = 1, add = TRUE)
plot(canals.x, col = "blue", lwd = 1, add = TRUE)
dev.off()
